# This script reproduces Tables 13 & 14 # 

rm(list=ls())
options(digits=4)
pkg <- c("PanelMatch", "ggplot2", "plm", 
         "matrixStats",
         "foreign", 
         "stringi", "stringdist",
         "DataCombine","lmtest", "multiwayvcov",
         "data.table",
         "dplyr", 
         "stargazer",
         "clubSandwich", "mediation",
         "doBy")
lapply(pkg, require, character.only = TRUE)

load("TG.RData")

## making a balance table
model1 <- lm(all_tour ~   
               age + 
               income + income_missing + male + 
               CCP + years_edu + internet + 
               as.factor(hukou) + 
               as.factor(occupation) + 
               as.factor(date) + 
               as.factor(city_county), 
             data = TG, 
             weights = postweight_agegender)
vcov1 <- clubSandwich::vcovCR(model1, type="CR1S", 
                              cluster = TG$SITECITY)
se1 <- coeftest(model1, vcov1)

model2 <- lm(all_tour_l31_180 ~   
               age + 
               income + income_missing + male + 
               CCP + years_edu + internet + 
               as.factor(hukou) + 
               as.factor(occupation) + 
               as.factor(date) + 

               as.factor(city_county), 
             data = TG, 
             weights = postweight_agegender)
vcov2 <- clubSandwich::vcovCR(model2, type="CR1S", 
                              cluster = TG$SITECITY)
se2 <- coeftest(model2, vcov2)

sink("output/Table13.tex")
stargazer(model1, model2,
          font.size = "large", digits = 4, 
          type = "latex",
          dep.var.labels = c("Current inspection", "inspection 30 to 180 days ago"),
          column.labels = c("Current inspection", "inspection 30 to 180 days ago"),
          covariate.labels = c("age", "income", "male", "CCP member", "years of education",
                               "internet usage"),
          dep.var.labels.include = F,
          style = "qje", no.space = FALSE, keep.stat = c("n", "rsq"),
          omit = c("city_county", "date", "income_missing",
                   "hukou", "occupation"
          ),
          omit.yes.no = c("No", "Yes"),
          add.lines = list(c("County FE?", "Yes", "Yes"),
                           c("Date FE?", 
                             "Yes", "Yes"),
                           c("Occupation FE?", 
                             "Yes", "Yes"), 
                           c("Household Registration FE?", 
                             "Yes", "Yes")),
          star.cutoffs = c(.10, .05, .01), 
          se = list(se1[,2], se2[,2]),
          p = list(se1[,4], se2[,4]), 
          header = F, notes = "robust standard errors clustered by prefecture in parentheses", 
          float = F)
sink()
###################3

# corruption satisfaction 
model1 <- lm(central_sat ~ all_tour + 
               age + 
               income + income_missing + male + 
               CCP + years_edu + internet + 
               as.factor(hukou) + 
               as.factor(occupation) + 
               as.factor(date) + 
               as.factor(city_county), 
             data = TG, 
             weights = postweight_agegender)
vcov1 <- clubSandwich::vcovCR(model1, type="CR1S", 
                              cluster = TG$SITECITY)
se1 <- coeftest(model1, vcov1)

model2 <- lm(central_sat ~ all_tour_l31_180 + 
               age + 
               income + income_missing + male + 
               CCP + years_edu + internet + 
               as.factor(hukou) + 
               as.factor(occupation) + 
              
               as.factor(date) + 
               as.factor(city_county), 
             data = TG, 
             weights = postweight_agegender)
vcov2 <- clubSandwich::vcovCR(model2, type="CR1S", 
                              cluster = TG$SITECITY)
se2 <- coeftest(model2, vcov2)

model3 <- lm(trust_in_ps ~ all_tour + 
               age + 
               income + income_missing + male + 
               CCP + years_edu + internet + 
               as.factor(hukou) + 
               as.factor(occupation) + 
               as.factor(date) + 
               as.factor(city_county), 
             data = TG, 
             weights = postweight_agegender)
vcov3 <- clubSandwich::vcovCR(model3, type="CR1S", 
                              cluster = TG$SITECITY)
se3 <- coeftest(model3, vcov3)

model4 <- lm(trust_in_ps ~ all_tour_l31_180 + 
               age + 
               income + income_missing + male + 
               CCP + years_edu + internet + 
               as.factor(hukou) + 
               as.factor(occupation) + 
               as.factor(date) + 
               as.factor(city_county), 
             data = TG, 
             weights = postweight_agegender)
vcov4 <- clubSandwich::vcovCR(model4, type="CR1S", 
                              cluster = TG$SITECITY)
se4 <- coeftest(model4, vcov4)


sink("output/Table14.tex")
stargazer(model1, model2, #model3, model4,
          font.size = "large", digits = 4, 
          type = "latex",
          dep.var.labels = c("Central Government Satisfaction", "Central Government Satisfaction"), 
          column.labels = c("Central Government Satisfaction", "Central Government Satisfaction"), 
        
          covariate.labels = c("current inspection", "inspection 30 to 180 days ago",
                               "age", "income", "male", "CCP member", "years of education",
                               "internet usage"),
          
          dep.var.labels.include = F,
          style = "qje", no.space = FALSE, keep.stat = c("n", "rsq"),
          omit = c("city_county", "date", "income_missing",
                   "hukou", "occupation", "Constant"
          ),
          omit.yes.no = c("No", "Yes"),
          add.lines = list(c("County FE?", "Yes", "Yes"), #"Yes", "Yes"),
                           c("Date FE?", 
                             "Yes", "Yes"), 
                           c("Occupation FE?", 
                             "Yes", "Yes"), #"Yes", "Yes"), 
                           c("Household Registration FE?", 
                             "Yes", "Yes") #"Yes", "Yes")
          ),
          star.cutoffs = c(.10, .05, .01), 
          se = list(se1[,2], se2[,2]), 
          p = list(se1[,4], se2[,4]), 
          header = F, notes = "robust standard errors clustered by prefecture in parentheses", 
          float = F)
sink()


